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ABSTRACT 

Lyman-limit absorption systems can play many important roles during and after cosmo- 
logical reionization. Unfortunately, due to the prohibitively large dynamic range required, it is 
impossible to self-consistently include these systems in cosmological simulations. Using fast 
and versatile semi-numeric simulations, we systematically explore the spatial distribution of 
absorption systems during and following reionization. We self-calibrate the resulting number 
of absorbers to the mean free path (mfp) of the ionizing ultraviolet background (UVB), and 
present results at a given mfp and neutral hydrogen fraction. We use a simple optical depth 
criterion to identify the locations of absorbers. Our approach is fairly robust to uncertain- 
ties such as missing subgrid structure. Unlike at lower redshifts where the UVB is relatively 
uniform, at higher redshifts the fluctuations in the UVB and the HII morphology of reion- 
ization can drive the large-scale distribution of absorption systems. Specifically, we find that 
absorbers are highly correlated with the density field on small scales, and then become anti- 
correlated with the UVB on large scales. After reionization, the large-scale power spectrum of 
the absorbers traces the UVB power spectrum, which can be predicted with a simple analytic 
extension of the halo-model. During reionization, absorbers tend to preferentially lie inside 
overdensities (i.e. filaments) of the recently-ionized intergalactic medium (IGM). Absorbers 
may also dominate the small-scale (k > 1 Mpc -1 ) 21-cm power during and after reioniza- 
tion. Conversely, they smooth the contrast on moderate scales. Once the HII regions grow to 
surpass the mfp, the absorbers add to the large-scale 21-cm power. Our results should prove 
useful in interpreting future observations of the reionization epoch. 

Key words: cosmology: theory - intergalactic medium - diffuse radiation - large scale struc- 
ture of universe - early Universe 



1 INTRODUCTION 

Recent years have been witness to impressive advances in com- 
put er modeling of the r eionization epoch (for a recent review, 
sec iTrac & Gnedini r2009). Current, state-of-the-art simulations of 
reionization are capable of including the interactions of dark mat- 
ter (DM), baryons, and an evolving ionizing radiation field on ~ 
100 /i _1 Mpc scales, while still resolving atomically-cooled ha- 
los at z < 10, which are thought to provide the bulk of the ion- 
izing photons in the adv anced stages of reionization l lTrac et all 
2008; Aubert & Tevssier 2010; for DM-only simulations see also 
McOuinn et alj|2007al ; ITrac & Cenll2007r). Such sim ulations have 
verified analytic models (e.g. iFurlanetto et al]|2004l) . which con- 
jectured that on large scales reionization proceeds in an "inside- 
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out" fashion, with overdense regions (which host the majority of 
sources) being ionized before underdense regions. 

Thus far, most of the advances in computer modeling have fo- 
cused on the sources of ionizing radiation. Current prescriptions for 
assigning ionizing luminosity to host halos and dealing with feed- 
back mechanisms are still rather crude. Justifiably, much research 
effort has gone (and continues to go) into exploration of these un- 
certainties (e.g. [McOuinn et alj|2007al ; iMesinger & Diikstr3l2008l ; 
ICroft & Al tav 2008). Nevertheless, we are quite confident in our 
ability to predict the distribution of DM halos that host the ionizing 
sources. 

Unfortunately, comparatively modest progress has been made 
in simulating the sinks of ionizing photons during reionization 
(with some notable recent exceptions; see below). These can po- 
tentially pla y important roles: as photon sinks, they can delay 
reionization dCiardi et ai]|2006l : [Alvarez & Abelll2010h : they can 
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affect the topology of reionization dChoudhurv et al.ll2009l) : they 
can dominate the ionizing photon mean free path (mfp) and so reg- 
ulate the ion izing ultraviolet background (UVB) during and afte r 
reionization jFurlanetto & Ohl2005l ; [Furlanetto & Mesingej2009h ; 
they can be responsib le for a non-negligible 21 -cm signal during 
and after reionization dChoudhurv et alj|2009l : IWvithe et all 2009a ; 
IWvithe & Loebll2009bl) . 

The difficulty lies in the dynamic range: small-scale structure 
would need to be resolved, along with the complicated feedback 
mechanisms regulating iQ Even at lower-redshi fts (z < 4) where 
there is a fair amount of observational da ta (e.g. IProchaskUl 1 9991 ; 
iPeroux et alj2003U ; |Prochaska et alj2009l) . consensus has failed to 
emerge on the nature of the dominant absorbers, the so-called Ly- 
man limit systems (LLSs), with studies associating them with low - 
mass dark matter halos (e.g. lKatz et al . 1996; Gardner etal.f2 001), 
galacti c or high-mass dark matter halos (e.g. iKohler & Gnedinl 
l2007ah or clumpy, so-called "cold-fl ows" of IGM gas streaming 
onto galaxies (e.g. iKeres et alj|2005l) . The situation is even more 
complicated at higher redshifts, where even more modest HI over- 
abundances can contribute significantly to the ionizing photon 
opacity. 

Instead of self-consistently simulating them, photon sinks in 
cosmological simulations are generally taken into account via so- 
called "clumping factors", which are a measure of t he mean re- 
combination rate inside sim ulation cells (though see Iciardi et al.l 
2006; iMcOuinn et al.ll20 07a). Clumping factors are often assumed 
to be spatially (and at times even temporally) uniform, though 
more s ophisticated techniqu es can use higher-r esolution simu- 
lations dKohleretal] |2007b) or analytic models IMcO uinn et al.l 
l2007ah when assigning clumping factors to cells. However, adding 
"by hand" these missing small-scale modes while preserving their 
large-scale distribution is very uncertain. Furthermore, because 
only the dumpiness of ionized gas is relevant for the recombina- 
tion rate, the clumping factor approach requires assumptions that 
( 1 ) the gas is already ionized and (2) that one knows a priori the dis- 
tri bution of the ionized and neutral regions (see e.g. the appendix 
in lFurlanetto & Ohll2005l ; llliev et al.l2005al) . 

Alternatively, some simulations have included discrete, self- 
shielded photon sinks in the f orm of collapsed, unresolved struc- 
tures v ia a subgrid approach dCiardi et al] 120061 ; IMcOuinn et al.l 
120074 These studies use the conditional excursion set formalism 
(e.g. iBond et al. I ll99ll : lLacev & Coldll993h to populate the simu- 
lation cells with unresolved collapsed objec ts, concluding that they 
can prolong the duration of reionization dCiardi et al . 2006) and 
have a sub-dominant role to sources in regulating the reionization 



1 Although current cosmological reionization simulations can resolve 
the Jeans length in the mean density, ionized IGM, most of the Universe 
was only exposed t o an ionizing background at redshifts around z ~ 10 
lLarson et alj .2010). Before then, the collapse of non-linear structures oc- 
curred in a neutral medium with a much smaller Jeans mass. As the small- 
scale structure of the baryons depen ds on the lo cal thermal history (not the 
instantaneous Jeans mass; lHui& Gnedinl 19971) . simulations would need to 
resolve the collapse of DM and baryons on extremely small scales (e.g. the 
Jeans length at z = 15 in even the mean density IGM is an impressively- 
small ^4 comoving pc), including self-consistent feedback mechanisms. 
The situation might not be quite so bleak, however, as some fiducial as- 
trophysical models predict an early epoch of X-ray heating (or some other 
form of heating), which would smooth fluctuations on T V \ T < fewxlO 3 
K scales, if the heating was early en ough and fairly homogeneous (e.g. 
iFurlane tto 2006a; Mesinger et al. 2010). Nevertheless, even if this dynamic 
range was achievable, a copious number of realizations would be needed to 
fully explore the dauntingly-large astrophysical parameter space. 



morphology dMcOuinn et alj|2007al) . Again, it is unclear if the dis- 
tributions of absorbers on large scales are well preserved in such 
subgrid models which only depend on the local cell's density. More 
fundamentally, it is unclear how to associate absorption systems 
with such collapsed "minihalos" (with virial temper atures < 10 4 
K) or with the halos hosting the ionizing sources (Wvi the et al.l 
l2009al) . 

In this paper, we take a more general approach. We do not 
attempt to model the progress of reionization and the accompany- 
ing role of absorption systems. Nor do we add missing small scale 
power in an effort to estimate the recombination rate inside each 
cell. Instead, we focus on predicting the large-scale distribution of 
abso rbers. 

iMiralda-Esctide et all (2000, MHR00) and subsequently 
Furlanetto & Oh ( 2005}) have assumed that absorption systems in- 
side the ionized IGM at z ~ 2—4 can be identified with a simple 
density threshold: regions of roughly the Jeans scale (computed 
at mean density and 10 4 K) whose density is above the thresh- 
old can be taken as fully neutral, while b elow the threshold, the 
IGM is fully ionized. In their appendix, Furla netto & Ob] (2005) 
have shown that this is similar to imposing a self-shielding cut- 
off, r > 1, at this critical density, given a homogeneous UVB. 
However, the ionizing background even in the cosmological HII 
regions becomes increasingly inhomogeneous at high er redshifts 
( Bolto n & Haehneltl2007l ; lMesinger & Furlanetto!: 2009), due to the 
increasing clustering of sources and the decreasing mean free path 
(mfp). 

Therefore in this paper, we assign self-shielded absorbers to 
cells with an optical depth criterion which takes into account the 
inhomogeneous UVB. We use semi-numericajf] simulations to ef- 
ficiently generate realizations of the density, halo, ionization (i.e. 
HII tomography), and ionizing flux fields. These are very fast 
compared to numerical simulations, allowing one to independently 
vary many parameters, and yet still produce HII morphologies 
which are in good agreement with radiative transfer simulations 
dZahn et alj|2007l; IMesinger & Furlanetto] 120071 ; IZahn et al]|2oTol : 
iMesinger et al] 20ld) . We stress that, lacking the dynamic range 
and radiative transfer, we do not presume to model the size and pho- 
toevaporation of absorbers, two important ingredients in answering 
how the absorbers affect the progress of reionization. Instead we 
present results at a given mean IGM neutral hydrogen fraction and 
mean absorber number density that we calibrate to the associated 
reduction in ionizing flux. 

Ou r approach is very sim ilar to the recent semi-numerical 
work o f lChoudhurv et al.l (2009). Motivated by the suggestion that 
reionization proceeds very gradually in a "photon-starved" regime 
(as su ggested by t he inte r pretation of the z > 5 Lya forest by 
Bolton & Haehnelt 200 7), Choud hurv et al] d2009l) have modeled 
the distribution of absorption systems and their impact on reion- 
ization morphology. To identify absorbers, they use an excursion- 
set approach comparing the local recombination rate to the time- 
integrated production rate of ionizing photons. Our approach dif- 
fers from theirs in several main points (see § [2] and § [3] for more 
details): (1) we estimate the reionization morphology and the in- 
stantaneous ionizing flux separately; (2) we use a self-shielded cri- 
terion instead of a recombination criterion to identify the location 
of absorbers, the latter being more sensitive to the missing subgrid 



2 By "semi-numerical" we mean using more approximate physics than nu- 
merical simulations, but capable of independently generating 3D realiza- 
tions. 
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small-scale power; (3) we use a self-calibrated, general approach to 
take into account (albeit crudely) missing/approximate physics. 

This paper is organized as follows. In §[2] we briefly summa- 
rize our semi-numerical modeling tools. In §[3] we present our pre- 
scription for identifying the locations of absorption systems. In §|4] 
we present our results, including statistical descriptions of the ab- 
sorber field after and during reionization. Finally in § 5, we sum- 
marize our conclusions. 

Unless stated otherwise, we quote all quantities in comov- 
ing units. We adopt the background cosmological parameters (Qa, 
f2 M , Q&, n, a 8 , H ) = (0.72, 0.28, 0.046, 0.96, 0.82, 70 km s -1 
Mpc -1 ), matching t he five-year results of the WMAP satellite 
dKomatsu et alj|2009h . 



z ~ 10. Note that homogeneous recombinations are implicitly in- 
cluded in the £ parameter. Also note that in this appro ach, unlike 
the recent modifications presented in IZahn et al.l 12010), the IGM 
is considered to be a two-phase medium, composed of either fully- 
ionized or fully-neutral cells. This is an excellent approximation on 
such small-scales, provi ded that stellar sp ectra dominate reioniza- 
tion (see the appendix in lZahn et alfeOlOl) . 

We ge nerate our ionizing flux fiel ds following the procedure 
described in Mesinger & Dijkstra (2008). We sum the contributions 
of halos with mass Mi at Xi for which the line of sight to the con- 
sidered cell at x only goes through ionized IGM: 



(1 + zf Mj X £ ion 

Z—i l-y — ir.\2 



-xi|/A 



4tt 



(1) 



2 SEMI-NUMERICAL SIMULATIONS 

We use the semi-numerical code Dexivff] to generate evolved 
density, halo, ionization and ionizing flux fields at z — 10. 
The halo, density and ion ization schemes are presented in 
iMesing er & Furlanetto ( 20 071, MF07), and the io nizing flux algo- 
rithm is presented in Mesinger & Diikstr3 d2008l) . We refer the in- 
terested reader to those papers for more details. Here we briefly 
outline these schemes. 

Our simulation box is L — 100 Mpc on a side, with the fi- 
nal density, ionization, and flux fields having grid cell sizes of 0.5 
Mpc. Halos are filtered out of the 1600 3 linear density field using 
excursion-set theory. Halo locations are then mapped to Eulerian 
coordi nates at a given redshift using first-order perturbation the- 
ory dZerDovichl 1970). Th is approach is s imilar in spirit to the 
"peak-patch" formalism of iBond & Mversl dl996l) . The displace- 
ment fields have an intermediate resolution of 800 3 , in order to 
most efficiently use the available RAM (velocity fields are corre- 
lated on larger scales than the density fields, thus they can afford 
worse resolution). The resulting halo fields have been shown to 
match both the mass function and statistical clustering properties 
of halos in N-body simulations, well past the linear regime (MF07, 
Mesinger et al., in preparation). 

The evolved (non-linear) density field is also computed in the 
same manner, by perturbing the 1600 3 Lagrangian density field 
with first-order perturbation theory. The resulting particle locations 
are then binned onto a Eulerian 200 3 grid. The statistical proper- 
ties of the density fields have been shown to match those from a 
hydrodynamic simulation remarkably well (Me singer et alj|2010h . 
Specifically, the density PDF agrees with the numerical simulation 
at the percent level for over a dex around the mean density at z > 7 
and roughly this cell size. Similarly, the density power spectra are 
extremely accurate at k < 5 Mpc -1 . We show slices through our 
density and halo fields in Fig.Q] 

To generate the ionizat ion field, we use excursion-set formal- 
ism ((Furlanetto et al. 2004) which compares the number of ion- 
izing photons produced in a region of a given scale to the num- 
ber of neutral hydrogen atoms inside that region. Specifically, we 
flag ionized cells in our box as those which meet the criterion 
/cou(x, z,R) ^ where £ is some efficiency parameter and 
/coll (x, z, R) is the collapse fraction smoothed around a cell at 
(x, z) on decreasing scales, R. The collapse fraction is computed 
using the resolved halo field, including halos down to a minimum 
mass of 10 8 M©, corresponding to the atomic cooling threshold at 
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Here, /(x, z) is the flux of ionizing photons at x (in s _1 pcm -2 ), 
and £ion(z) is the the rate at which ionizing photons are released 
into the IGM by a dark matter halo per unit mass. We assume 
a fiducial value of e- lon {z) — 3.8 x 10 5S fif,/ Tly!_/tu{z) photons 
Mq 1 s -1 which provides a good fit to the obser ved z ~ 6 lumi- 
nosity functions of Lyq emit t ing galaxies (LAEs) dKashikawa et al.l 



20061-IShimasaku et al] 2006; Dij kstra et alj2007l ; IStark et al.l2007l; 
McOuinn et alj l2007bf). and the z = 6 Lyman Break galaxies 
(LBGs) (e.g. lBouwens et"aT1l2006l) . However this fiducial normal- 
ization can be subsumed in our optical depth normalization (see 
§ [3j and is not relevant when presenting results at a fixed number 
density of absorbers. What we are really after is the spatial distribu- 
tion of flux, and its impact on the distribution of absorbers. Hence, 
where applicable we present most results in terms of the normalized 
flux, 



/n(x,z) = f(x,z)/f(z) 



(2) 



where / is the mean value of the flux field for xm ~ 0. We chose 
to normalize by the flux post-reionization to facilitate direct com- 
parisons at the various reionization stages. 

The free parameter A in eq. |Q3 is the assumed mean free path 
of the UV photons in the ionized IGM. We use the mean free path 
to self-calibrate the number of absorbers with the flux fields, as de- 
scribed in § [3] There is considerable disagreement among present 
observational estimates of the mfp at the highest redshifts acces- 
sible to such observations, z ~ 4 (Storrie-Lombardi et al.ll 19941; 
I Stengler-Larrea et al.l [l995t IPeroux et al.l l2003bl Prochaska et al.l 



I2009L c.f. their Fig. 14. ISongaila & Cowid 120101) . Furthermore^ 
there is no reason to expect that the empirically-derived redshift 
evolution of A in these studies should extend to higher redshifts. 
Therefore, we try to be general in our conclusions, and show re- 
sults from two fiducial values of A = 10 and 20 Mpc. 

For our A = 10 (20) Mpc cases we obtain / ~ 2.5 (4.3) x 10 4 
ionizing photons s -1 pcm -2 , which corresponds to the mean ion- 
ization rate f ~ 2 (3.3) x 10 -13 s -1 . Note that in addition 
to matching the LAE and LBG luminosity functions as men- 
tioned above, our procedure results in values of T which are 
consistent with recent analysis of the z ~ 5-6 Lyman-a for- 
est, provided F does not evol ve much from z ~ 6 to 2 ~ 10 
(e.g. iBolton & Haehnell l2007l ; note that F already was found not 
to evolve much from z ~ 3 to z ~ 6). For more detailed 
analyses of the ionizing fl ux distributions, we refer readers to 
IMesinger & Diikstral d20o3) . 

Finally, we stress that ou r flux fi elds, unlike those implicitly 
computed in IChoudhurv etafl (2009). take into account the com- 
plex morphology of the HII regions. HII bubbles ar e not very spher- 
ical in the advanced stages of reionization (e.g. iMcOuinn et al.l 
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Figure 1. Slices through the matter overdensity field, <5(x, z) = p/p — 1 (left panel) and the corresponding halo field (right panel). Both slices are 100 Mpc 
on a side and 0.5 Mpc deep. 



l2007al : lTrac & CerfeOOl MF07). Therefore, many sources within 
the mfp will still not contribute to the local ionization rate, since 
they do not have a "clear line of sight" to the cell. Instead, these 
photons will go into expanding the HII regions and altering their 
morphological structure. 



3 IDENTIFYING SELF-SHIELDED ABSORPTION 
SYSTEMS 

Using the modeling tools summarized in the previous section, we 
compute the density, halo, HII morphology, and ionizing flux fields. 
The next step therefore is to identify the likely locations of self- 
shielded absorption systems. As mentioned in the introduction, we 
do this with an optical depth criterion, which takes into account the 
inhomogeneity of both the density field and the UVB. Simulation 
cells with 

t(x,z)>1 (3) 

are estimated to contain enough neutral gas to host self-shielded 
absorption systems. Note that since we present results at a fixed 
number density of absorbers, the precise value of the RHS of eq. l[3) 
is irrelevant (see below). The optical depth across a cell is estimated 
according to 

t(x, z) = A XmriHcradx , (4) 

where da; is our cell length of 0.5/(1 + z) pMpc, and <jh ~ 
6 x 10~ 18 pern 2 is the photoioni zation cross section of the hydro- 
gen dOsterbrock & Ferlanj|2006l) , n H (x,z) is the cell's (proper) 
number density of hydrogen, xm is the local neutral fraction com- 
puted assuming ionization equilibrium, 

nHZHiT = n|(l — xm) 2 ct B , (5) 



wh ere qb is the case-B recom bination coefficient at T = 10 4 
K dOsterbrock & Ferlandl2006h , and F(x, z) = /(x, z) <th is the 
local ionization rate computed using eq. ([TJ. 

The parameter, A, in eq. I0 is used to roughly (and homoge- 
neously) incorporate missing/uncertain physical components of our 
approach (e.g. the absorber cross-section; see below). We use it to 
self-calibrate the resulting number density of absorbers according 
to 

A= [V2n ss {A)dx 2 Y 1 , (6) 

where nss {A) — Nss (A) / L 3 is the number density of the self- 
shielded regions and Nss {A) is the number of self-shielded pixels 
found through the self-shielding criterion. We self-calibrate our re- 
sults by choosing A such that the RHS of eq. {6]l matches the mfrjf] 
used to calculate the flux fields, i.e. eq. (TJ. For better statistics, this 
calibration is done in a fully-ionized universe; therefore the free- 
parameter A is the same for all values of xhi and only depends on 
the mean separation of absorber jfl 

Our cell size roughly matches the Jeans mass in the ionized 
IGM, thus crudely mimicking the effects of Jeans smoothing. How- 
ever, as noted above, small-scale fluctuations are not set by the in- 
stantaneous Jeans mass, but instead depend on the local thermal 



4 The Case A value is more appropriate if the ionizing photons are ab- 
sorbed inside dense, neutral systems (e.g. MHR00); however, at higher red- 
shifts, an increasin g fraction of photons ge t absorbed in moderate optical 
depth systems (e.g. Furlanetto & Oh 2006c). In any case, the choice of re- 
combination coefficients is not important, as it gets subsumed in the "A" 
normalization parameter. 

5 Note that eq. (6) directly normalizes the resulting number density of ab- 
sorbers, not the mfp. If absorbers are inhomogeneously distributed, the ac- 
tual mfp can be different than predicted by the above equation. 

B Since A is related to the intrinsic self-shielding properties of our sim- 
ulation cells, it should be much more sensitive to the mfp than to xjji- 
Furthermore the above calibration only makes sense if the attenuation of 
the ionizing background is indeed dominated by the number density of ab- 
sorbers; i.e. when the ionizing photons are no longer absorbed by the HII 
bubble walls. This is only the case in the advanced s tages of reionization, 
when the HII bubbles grow to be la rger than the mfp tMes inger & Diikstral 
l2008t l iFurlanetto & Mesinger 2009). In the early stages of reionization be- 
fore this transition, the ionizing background is instead dominated by the HII 
morphology and is therefore insensitive to the mfp in the HII regions. 
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history and their state of evolution when first exposed to a UVB. 
Therefore it is very likely that, like all cosmological reionization 
simulations, we are missing subgrid structure. Processes such as 
the subgrid recombination rate and the photoevaporation rates of 
absorbers are indeed very s ensitive to the detailed profiles of the 
ionized and ne utral gas (e.g [iiiev et al.ll2005b| ; ICiardi et al]|2006l ; 
McOu rnn et al.l l2007al : IChoudhurv et al.l |2009|) . However, in this 
paper we are merely concerned with identifying the locations of 
self-shielded systems. For this task, the simple criterion of "there is 
enough neutral gas in the cell" of eq. ([3} is likely a decent indicator 
that there are o ne or more absorption systems inside the cell (c.f. 
the appendix of Furla netto & Ohl2 005), Furthermore, any homoge- 
neous dependence of our results on the missing small-scale power 
(such as the absorber cross-section) is incorporated via our calibra- 
tion parameter A. We obtain A ~ 3 (4) for our A = 10 (20) Mpc 
models. 

Unfortunately, there is an inconsistency in the above self- 
calibration. As we shall see below, the distribution of absorption 
systems at high redshifts is likely very inhomogeneous. This makes 
the smooth flux attenuation term, exp[— r/A], only a rough approx- 
imation. Such a smooth attenuation is more naturally interpreted 
as being due to the integrated column density of low-r systems 
(i.e. the Lya forest). If these more common systems did in fact 
dominate the absorption of ionizing photons at high redshifts, then 
our assumption of a global mfp is quite reasonable. However, if 
the flux was dominated by such low-r systems, then our calibra- 
tion in eq. l[6j would be unnecessary, since the flux is not sensitive 
to self-shielded absorbers. Then our results would best be inter- 
preted at a fixed number density of absorbers, with the absorber 
abundance unrelated to the attenuation of the UVB. Interestingly, 
iFurlanetto & O3{2005) find that at z > 6, the contribution of self- 
shielded LLSs and the Lya forest to the UV optical depth is roughly 
comparable, though this estimate is model-dependent and only ap- 
plied to the post-reionization regime. Indeed, since our values of 
A are somewhat higher than unity, it is possible that some of the 
absorption systems we flag are in fact lower-r cells more appropri- 
ate for the Lya forest (see eq.@. An improvement to our approach 
would be to assign parameterized cross-sections to each absorber 
and iteratively compute the UVB (which could have additional 
smooth attenuation from low-r systems) and absorber fields until 
the two converge. In this way the UVB and absorber distributions 
would be self-consistent. However this is both time-consuming and 
of dubious benefit, given our present ignorance of absorber prop- 
erties at high-redshift. Given the above uncertainties, one should 
keep in mind that our quantitative estimates below are less robust 
than our qualitative ones. 

Although we allow absorbers to be hosted by the same cells 
hosting ionizing sources, we do not explicitly model the HI con- 
tent of the interstellar medium (ISM). Resolution requirements and 
detailed feedback processes make it very difficult to simulate a sta- 
tistical sample of such absorbers. At lower redshifts, high opac- 
ity systems such as Damped Lyman alpha systems (D LAs) are in- 
deed associated with galactic mass dark matter halos I Wol fe et"al] 
2005). On the other hand, low opacity systems such as the Ly- 
man a fo rest are assoc iated with the filamentary structure of the 
web (e.g. ISchavduOOll) . At z ~ 5, the to tal UV opacity is dom- 
inated by the inte rmediate opacity LLSs (Miralda-Escude 2003; 
IPeroux et alj2005l) , with the nature of these systems being the least 

understood. As already mentioned, the IGM can self-shield much 

J 1 1 n 

easier at higher redshifts ( IFurlanetto & Oh 2005). Therefore it is 

likely that at the epochs we consider, the total UV opacity is in fact 



dominated by absorbers residing in the IGM, which are simpler to 
model. 

Finally, we point out that our method, unlike the one in 
Choud hurv etafl (2009), computes the reionization morphology 
separately from the absorber locations. The overall reionization 
morphology depends on the complete ionization history, whereas 
the absorber distributions are more sensitive to the instantaneous 
ionizing flux. Ind eed, using large - scale, DM-only simulations with 
radiative transfer McOuinn et alj d2007al) find that absorption sys- 
tems have a small impact on reionization morphology, even with 
extreme assumptions, such as a small mf p and no photoevapora - 
tion (see the bottom panel of Fig. 10 in iMcOuinn et"all l2007a). 
Similarly, shadowing effects from the absorbers should be small 
since ionizing sources are comparably very abundant with many 
sources from many directions contributing to the local ionizing 
background. Nevertheless, we stress that treating absorption sys- 
tems separately from the HII morphology is an assumption of our 
method, not a result. 



4 RESULTS 

4.1 Absorber Distributions After Reionization 

Here we briefly discuss our results post-overlap, i.e. when thi = 
o3. This is the simplest case we consider, since the HII morphology 
has no impact on the results. Therefore it is a good place to start and 
build some physical intuition, before proceeding to the reionization 
epoch. 

In Figures [2] and [3] we show slices through the flux (top pan- 
els) and ionization (bottom panels) fields, generated assuming A = 
10 Mpc (Fig.|2j and A = 20 Mpc (Fig.[3}. The right-most columns 
correspond to the post-reionization regime. 

We see that even post-reionization, the UVB is spatially in- 
homogeneous, with the inhomogeneity increasing with decreas- 
ing mfp. Even on large, ~10 Mpc scales, spatial fluctuations in 
the UVB can be a factor of a few around the mean; on smaller 
scales, fluctuations can even be a factor of ten around t he mean 
dMesinger & Diikstrall2008l ; iMesinger & Furlanettdl2009l) . This is 
in sharp contract to the moderate and low redshift regime (z < 
4), where spatial fluctuations contribute o nly a few percent of 
the mean value of the UVB Jziiol 1 19921; iFardal & Shul 
Gnedin & Hamiltonl |2002| ; iMeiksin & White! |2004| ; ICrofl 
Bolton et al. 200j|- 



1993 



2004 



These fluctuations in the UVB indeed have a strong impact 
on the distribution of absorbers. Looking at the corresponding dis- 
tributions of absorbers (shown as green dots in the bottom pan- 
els of Figures [2] and [3}, one can see that they preferentially lie 
in regions with a weaker UVB, far from large clusters of ionizing 
sources (shown as black dots). This effect is even more pronounced 
when the mfp is larger. As a result, the absorbers themselves are 
spatially inhomogeneous. If this absorber inhomogeneity was self- 
consistently included in our UVB calculation, it would further in- 
crease its fluctuations. 

On the other hand, some absorbers are seen to lie close 



7 We present results in this work at fixed Shi. which is computed without 
the contribution from the HI inside the absorption systems. This is the most 
natural way to highlight the additional impact of absorption systems, and 
it also avoids having to assume how much of the absorber's gas is actually 
neutral. 
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Figure 2. Upper panels: the normalized flux field calculated assuming A = 10 Mpc, for xhi = 0.72, 0.45, 0.18, (left to right). Lower panels: the corresponding 
ionization fields. Neutral regions are shown in red, while ionized regions are shown in white. Halo locations are marked with black points, while absorber 
locations are marked with green points. All slices are 100 Mpc on a side and 0.5 Mpc deep and correspond to the same region shown in Fig.[T] 




Figure 3. Same as Fig. [2] but assuming A = 20 Mpc. 



to ionizing sources. These absorbers likely lie in overdense re- 
gions, where the high gas density manages to shield them from 
the increased ionizing radiation. Note that the optical depth scales 
roughly as r oc (8 + l) 2 / -1 , 

To quantify these trends in Fig. [4] we show the cross- 
correlation function, £(-R), of the binary absorber field with the 



density (left panel) and the flux (right panel) fields. This figure 
confirms the qualitative intuition from above. At small scales, ab- 
sorbers are highly correlated with overdensities. This correlation 
increases as the background becomes more uniform and there are 
fewer absorbers. However, on large scales, absorbers become anti- 
correlated with the density field instead lying in large-scale under- 
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Figure 4. The cross-correlation between the self-shielded regions and 8 (left panel) or fjy (right panel), at Shi = 0. The blue solid and the red dashed curves 
correspond to A = 10 and 20 Mpc, respectively. 



dense regions where the UVB is weak. This effect is even more 
prominent in the strong anti-correlation of the absorbers and the 
UVB on large-scales. Note that our boundary conditions which 
suppress modes larger than the box size bring the cross-correlation 
back to zero as one approaches the box size. 

Therefore scales are very important in understanding the dis- 
tribution of absorbers. On small scales, the density dependence of 
the optical depth dominates, and the density field determines the 
absorber locations. On the other hand, on large scales it is the ion- 
izing flux which governs the locations of absorbers. This is ex- 
pected since the UVB is correlated on larger scales than the density 
field. The scale at which this transition happens and the strength 
of the (anti-)correlation depend on the mfp. Thus, our results argue 
against the common practice of adding substructure to ionization 
cells based on just the local density; long- wavelength modes should 
also be considered. In particular, the sources are much more biased 
than the absorbers, so fluctuations in the ionizing background - 
which extend to large scales around a clump of sources - dominate 
the variations in the absorber population. 



4.2 Absorber Distributions During Reionization 

Now that we have a better understanding of the scale-dependent 
absorber distributions post-reionization, we explore how these are 
affected by the HII morphological structure during reionization. 
Reionization morphology further modulates the ionization field, as 
can be seen in the top panels of Figures [2] and [3] From the bottom 
panels of these figures, we confirm that the absorbers preferentially 
cluster in regions of weak flux. During reionization, this translates 
to the edges of HII regions.When the process is viewed at a fixed A, 
the absorbers essentially "move-in" to the newly-ionized regions. 
However, inside those large-scale regions, absorbers will preferen- 
tially be located in overdensities, i.e. filaments (see also Fig.[8j(. 

These trends can also be seen more quantitatively in Fig. [5] 
where we plot the cross-correlation of the halo and absorber fields 
for Shi = 0.72, 0.45, 0.18 and 0, left to right. The halo field is 
computed on the same resolution, 200 3 grid, with cells containing 
halos flagged as T, and those not containing halos flagged as '0'. 
We see that early in reionization, when the HII regions are smaller 
than the mfp, there is only a positive correlation between the halos 
and absorbers. This is understandable since by definition absorbers 
must reside inside HII regions, and on these relatively small scales, 
they both correlate to the density field. However, as the HII regions 



grow, the large-scale fluctuations in the UVB become important in 
determining the absorber locations. Absorbers show an increasing 
preference for the newly ionized regions distant from the ionizing 
sources where the UVB is weak. This manifests itself through the 
increasingly negative large-scale cross-correlation as reionization 
progresses. 

In Fig. [6] we plot the distribution of the absorption systems 
in the (5,/jv) plane. The black, blue and red surfaces enclose 99%, 
95% and 68% of the absorber population, respectively. The upper 
(lower) panels correspond to A = 10 (20) Mpc. This figure clearly 
shows the degeneracy between the density and ionizing flux, as the 
contours stretch from the lower left to the upper right, following 
r oc (S + l) 2 // isocontours. The isocontours are flatter for higher 
values of the mfp due to our choice of normalization of the ion- 
izing flux, which means that the same /jv corresponds to a higher 
unnormalized value of the UVB for A = 20 Mpc than for A = 10 
Mpdl 

Fig. |6]is further confirmation of our narrative above. Early in 
reionization, the highly biased nature of the ionizing sources im- 
plies that ionized regions generally trace large- scale overdensities , 
i.e. reionization is inside-out on large scales dTrac & Cenl 120071 ; 
iMcOuinn et"al]|2007al ; IZahn et alJboiCh . Within these large-scale 
regions, it is the density fielcQ which determines whether a cell 
will host absorption systems; hence the contours in the left pan- 
els of Fig. [6] are fairly steep showing insensitivity to the UVB. As 
the HII regions grow to include large-scale voids, it is the UVB 
which increasingly determines whether a cell will host absorption 



8 Note that our normalized flux, /jy from eq. (2)> is normalized to the 
mean value post reionization, for each mfp. Therefore the same /jv values 
correspond to the same unnormalized ionizing flux for each Shi. but are 
different for different values of A. The mean value of the UVB is of course 
higher for higher values of the mfp. 

9 Note that the absorber densities in the figure are much more modest (of 
order unity) than in the post-reionization models of MHR00. MHR00 stud- 
ied self-shielded systems at lower redshifts, z = 2-4. At redshift z = 10, 
the mean density of the Universe is an order of magnitude greater, and the 
IGM can self-shield at much smaller overdensities since r oc A 2 . Also, 
we study much smaller mean free paths, which effectively means that the 
UVBs that we consider are much weaker, and absorbers are more com- 
mon. Finally, our semi-numerical simulations are at a lower resolution, so 
we have less small-scale power than MHR00. We remind the reader that 
our procedure flags the locations of the absorbers but does not necessarily 
require that the absorbers themselves fill the entire cell. 
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Figure 5. The cross-correlation between the absorption systems and the halo field at z = 10, for all the neutral fractions considered in this work. The blue 
solid and the red dashed curves refer to A = 10 and 20 Mpc, respectively. 
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Figure 6. The distribution of the absorption systems in the (<5,/jv) plane. The black, blue and red surfaces enclose 99%, 95% and 68% of the absorber 
population, respectively. The upper and the lower panels correspond to A = 10 and 20 Mpc, respectively. 



systems; hence the contours flatten as reionization progresses. As a 
consequence, the absorbers "move-into" the filaments of the newly- 
ionized large-scale underdense regions. This can be seen as the shift 
in the contours towards higher density cells. Note that cells occu- 
pying the lower right quadrant of the figure still host absorbers, but 
that region of parameter space is too rare to be enclosed by the 99% 
contour. 



Although a direct comparison with Choudh urv et alj 12009) is 
difficult since they do not explicitly distinguish between the neu- 
tral IGM outside HII regions and the absorbers inside HII regions, 
they find a similar qualitative trend of the neutral IGM remain- 
ing in increasingly overdense regions as reionization progresses. 
However, we do not find evidence for their dual peak in the den- 
sity PDF of the neutral IGM at J > 0. Althoug h the less-dense 
neutral regions found by Chou dhurv etafl (2009) could represent 
the neutral IGM outside HII regions, they still obtain significant 
neutral cells with 8 > 0, which is not consistent with the inside- 
ou t nature of reionization on large-scales. The parameter studies 
bv lMcOuinn et al . (2007a) suggest that absorbers have a small im- 



pact on reionization morphology, and that the inside-out nature of 
reionization on large-scales is a robust result. As already explained, 
this was one of the reasons we separately compute the reionization 
m orphology and ab s orber locations, unlike the analytic procedure 
of IChoudhurv et al.l (2009) which combines the two calculations. 
We stress that when the UVB is inhomogeneous, the distribution of 
absorbers is best represented in the (5,/jv) plane, and not marginal- 
ized over just the density. 

The above-mentioned trends can also be seen Table [T] where 
we show the fraction of cells containing absorbers, normalized by 
the HII filling factor (x^), and unnormalized (s a ). Indeed the 
abundance of absorption systems inside the ionized IGM, shown by 
a; a N is much higher in the early stages of reionization. Given that 
reionization is inside-out on large scales, this lends further confir- 
mation that the density field dominates absorber locations early in 
reionization. 

In Fig. [7] we show the power spectra of our various fields. The 
dimensionless power spectrum of field i is taken to be A 2 (k) = 
k 3 /(2n 2 V) (\Si(k)\ 2 }k, where $i(x) is the dimensionless fluctu- 
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^'HI 


A/Mpc 




■'a 


0.72 


10 
20 


0.22 
0.27 


0.06 
0.08 


0.45 


10 
20 


0.14 
0.14 


0.08 
0.08 


0.18 


10 
20 


0.09 
0.05 


0.08 
0.04 





10 

20 


0.03 
0.02 


0.03 
0.02 



Table 1. The fraction of simulation cells containing absorbers, normalized by the HII filling factor (:EaNX an d not normalized by the HII filling factor (x a ). 



ation of field i. The top panels correspond to the density and the 
UVB. The bottom panels correspond to the absorber power spectra 
at x m = 0.72, 0.45, 0.18 and 0, left to right, with A = 10 (20) Mpc 
shown with the blue solid (red dashed) curve. 

Since absorbers strongly correlate with the density field on 
small scales, we would expect their small-scale power spectrum to 
roughly follow the density power spectrum, with some bias. On 
large scales however, the absorbers anti-correlate with the flux (and 
density) fields. Therefore, their large-scale power should follow the 
UVB power spectrum, once HII bubbles have grown large enough 
to allow the absorbers to populate these modes. We see that these 
trends are indeed present in the power spectra in Fig. [JJ There is a 
dip in the absorber power spectrum between these two limits where 
the absorbers transition from being correlated to the density field 
on small scales to being anti-correlated on larger scales. 

Differences in the absorber power spectra of our two mfp mod- 
els only start emerging when the HII bubbles grow large enough to 
surpass the mfp (the mean HII bubble scalqlj is shown with the 
black vertical line in Fig. [7j. The A = 20 Mpc curves are in gen- 
eral higher during reionization, since there are fewer absorbers, and 
so the fractional offsets shown as the power spectra in Fig. [JJ are 
higher. When shown in dimensional units (scaled by the squared 
mean signal), a smaller mfp does have more power, as we shall see 
in 5l4~3l 

Unfortunately, our simulation box is only 100 Mpc on a side, 
which is only a few times the mfp. Hence, it is difficult to fully iso- 
late the impact of A from the absorber power spectra. To gain fur- 
the r insight, we turn to an analytic model for the UVB power spec- 
tra. lMe singer & Furlanettol {2009) have shown that the UVB power 
spectra can be predicted using an analogy of the "halo model", i.e. 
as the sum of two terms, one arising from correlations within a sin- 
gle source's 1/r 2 envelope, and o ne arising from correlations be- 
tween the locations of two sources dScherrer & Bertschingej[T99ll ; 
ICoorav & S heth 2002). The large-scale power corresponds to the 
two-source term, which traces the linear matter power spectrum on 
large scales times the bias, before being truncated at k > 1/A. We 
plot this two source term as the magenta curves in the flux panel 
of Fig. [JJ The analytic model confirms that the A = 10 and 20 
Mpc curves indeed do converge on large scales. The peak and drop 
at k > 1/A (demarcated by the blue and red vertical lines) is also 
seen in the absorber power spectra at xm = 0. The detection of 



this feature would be a novel measurement of the ionizing pho- 
ton mean free path, valid even if the mfp was not dominated by 
the self-shielded absorbers themselves. We defer a detailed ana- 
lytic treatment of the absorber power spectra, including cross- and 
higher-order terms, to future work. 



4.3 Impact of the Absorbers on the 21-cm Signal 

Some of the most important and timely observations of the epoch 
of reionization will come in the form of the redshifted 21-cm 
line of neutral hydrogen. Several interferometer telescopes are 
scheduled to become active in the near future, including the 
Mileura Wide Field Array (MW A; iBowman et al.1 [^OO SPl the 
Low Frequency Array (LOFAR; lHarker et alj I201CIPI, th e Gi- 
ant Metrewave Radio Telescope (GMRT; |Pen et alj|2003) . the 
Precision Array to Probe the Epoch of Reionization (PAPER; 
IParsons & PAPER Tearj2 009). and eventually the Square Kilome- 
ter Array (SKAQ. The neutral hydrogen of the absorbers could 
have a significant impa ct on the interpre t ation of the reionization 
signal, as suggested bv lChoudhurvet alJ (2009). Here we investi- 
gate this signal in our models. Note again that we do not explicitly 
focus o n the HI content of the ionizing galaxies with T v i r > 10 4 
K, as in lWvithe et all d2009j> . 

The offset of the 21-cm brightness temperature from the CMB 
temperature alon g a line of sight (LOS) a t observed frequency v can 
be written as (c.f. lFurlanetto et alj|2006bl) : 



ST b (x,z) 
mK 



9(l + z) 1/2 x m (x,z){l + S nl (x,z)] , (7) 



where we have chosen to ignore peculiar velocity effects, shown 
to have a rel atively minor impact during t he advanced stages o f 
reionization I McOui nn et al ] |2006L MF07, lMesingeretal.ll2010h . 
Equation {7} further assumes that the spin temperature of the IGM 
has been heated well past the CMB temperature, which seems 
likely durin g the bulk of reionization for fiducial astrophysical 
models (e.g IPritchard & Furlanet to 2006; Furla netto et al.|[2 006b; 
Chen & Miralda-Escudd 120081 ; ISantos et alJ 120081 : iBaeketaU 
20091; iMesinger et alfcOld) . 

Slices through our 21-cm brightness temperature boxes are 
shown in Fig.U for x m = 0.72, 0.45, 0.18 and (left to right) and 
A = 10 (upper panels) and 20 Mpc (lower panels). Black regions 



10 The mean bubble size is calculated by randomly choosing an ionized 
cell and recording the distance from that cell to the edge of the HII region, 
along a randomly chosen direction. These distances are then averaged to- 
gether to obtain the mean bubble size. 



11 http://web.haystack.mit.edu/arrays/MWA/ 

12 http://www.lofar.org 

13 



http://www.skatelescope.org/ 
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Figure 7. Top panels: the power spectra of the matter fluctuations <5 (left) and the UVB (right). The black, blue, green and red curves refer to 5'hi = 
0, 0. 18, 0.45 and 0.72, respectively, for A = 1 (solid lines) and 20 Mpc (dashed lines). The magenta lines in the right panel correspond to the im = 
analytic model of Mesinger & Furlanetto (2009) in the large-scale limit; see text for details. Lower panels: the power spectra of the absorbers, for A = 10 
(blue solid lines) and 20 (red dashed lines) Mpc. The colored vertical lines correspond to A; = 1/A, while the black dot-dashed vertical line represents the 
mean size of the HII regions (computed according to MF07). 




Figure 8. Slices through our <5Tj, boxes, assuming A = 10 (upper panels) and 20 Mpc (lower panels), and xhi = 0.72, 0.45, 0.18 and (left to right). All 
slices are 100 Mpc on a side and 0.5 Mpc deep and correspond to the same region shown in Fig.[T] 



correspond to the HII bubbles, with the absorbers seen to preferen- 
tially populate the outskirts of these regions. The filamentary struc- 
ture of the neutral IGM can be seen outside the HII regions, with 
the strongest (reddest) signal coming from overdensities. Looking 
closely at Fig.[8]one can note that indeed as the HII regions expand 
into the neutral IGM, the absorbers preferentially correspond to the 
overdensities in the newly-ionized medium. 

The corresponding dimensional 3D power spectra, 
5T b 2 A 2 21 {k), are shown in Fig. with A = 10 (20) Mpc 
corresponding to the blue solid (red dashed) curves. The black 



dot-dashed curve corresponds to the 21 -cm power neglecting 
absorption svstemj^l 



1 As mentioned previously, Shi is computed without the contribution 
from the HI inside the absorption systems. The "total" neutral fraction of 
the Universe would be xj^j = Sjji + x ^\^ wnere x hi is the fraction of 
the IGM gas contained as HI inside absorption systems. In our fiducial case 
where all of the absorber's gas is assumed to be neutral (an assumption only 
relevant for the quantitative predictions in this section), x^f ~ which 
can be read from Table 1 . 
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Figure 9. The dimensional power spectrum of 21-cm fluctuations, estimated assuming A = 10 (blue solid lines) and 20 Mpc (red dashed lines), for xhi = 
0.72, 0.45, 0.18 and (left to right). The black dot-dashed curve corresponds to the 21-cm power neglecting absorption systems. 



Neglecting absorbers, this signal progresses along the familiar 
story lines. Before reionization, the 21-cm power spectrum is dom- 
inated by the matter power spectrum. Reionization carves out dark 
regions in the 21-cm signal, which initially (Shi > 0.9) steepen the 
power spectrum as the first regions to be ionized ar e the small-scale 
overde nsities, equilibrating power on large scales dMcOuinn et al.l 
l2007al) . As the HII regions grow, the power spectrum is governed 
by HII morphology, with a "shoulder" feature corresponding to the 
bubble scale pro gressing from smal l to large scales, and stalling at 
the mfp scale iMesinger et alj[201ol) . 

Including absorbers complicates this story. Most notably, they 
drive small-scale power. Therefore the 21-cm power spectrum in- 
creases at small scales (k > 1 Mpc~*) where the HII bubbles 
would have otherwise damped the fluctuations. On moderate scales 
(0.5 Mpc -1 < k < 1 Mpc -1 ), the absorbers actually smooth the 
fluctuations somewhat by decreasing the contrast between the ion- 
ized and neutral regions. Once the HII regions grow to be larger 
than the mfp, the absorber distribution (sourced by the UVB) adds 
to the large-scale power. 

Again we note that the power spectrum is not sensitive to the 
mfp until the later stages of reionization. The modulation of the 
UVB by a smaller mfp and the increased number density of ab- 
sorbers result in a stronger signal of the A = 10 Mpc model than 
the 20 Mpc model. 

As noted by Choud hurv et~all d2009h . the IGM absorbers 
can c ontribute to a sizable 21-cm signal (though see llidz et al.l 
2008) during and po st-reionization (more so than galactic HI; e.g. 



Wvit he et al.ll2009ah. Overall, we find a smaller 21-cm signal than 
Choud hurv et alj ( l2009h . who also do not find the non-monotonic 



behavior of the power spectrum noted above. If indeed the 21-cm 
signal from absorbers is strong even post-reionization, this could 
bias the calibration of upcoming interferometers. 

We remind the reader that our procedure assigns all of the gas 
inside the absorber's simulation cell as being fully neutral. This 
need not be the case, as mentioned previousl y. Indeed, low redshif t 
LLSs are found to be mostly ionized (e.g. IPeroux et alj|2003bl) . 
although absorption systems during the reionization epoch likely 
have different properties. Nevertheless, readers should note that the 
absorbers could contribute less to the 21-cm signal than shown 
henB In fact, if we allow the absorber host cells to be partially 
ionized, a residual HI fraction of only 10 -3 - 10 -2 is sufficient 



15 In the previous sections we focus on the spatial distribution of ab- 
sorbers, which should be robust to these uncertainties. Furthermore, the di- 
mensionless power spectra in Fig.|7]is less sensitive to subgrid structure, 
since the fractional variation, S = pip — 1 is unchanged for a constant 
absorber ionization fraction. 



to achieve a lyman limit optical depth of r ~ 1. This would de- 
crease the absorber contribution to the 21-cm power spectrum by 
roughly the same amount, making it negligible. Partially ionized 
high-z absorbers might be more common if, for example, higher 
energy photons can efficiently penetrate the absorbers. Regardless 
of what exactly is the neutral fraction of the absorber cells, the qual- 
itative trends in the 21cm power spectrum discussed in this section 
should remain unchanged. 



5 CONCLUSIONS 

Lyman-limit absorption systems can play many important roles 
during and after cosmological reionization: as photon sinks, they 
can delay reionization; they can regulate the UVB; they can affect 
the morphology of reionization; they can be responsible for a non- 
negligible 21-cm signal. Unfortunately, due to the prohibitively 
large dynamic range required, it is impossible to self-consistently 
include these systems in cosmological simulations. 

In this paper, we focus on predicting the spatial distributions 
of absorbers during and following reionization. We use the efficient 
semi-numerical simulation code, DexM, to generate density, halo, 
ionization, and UVB fields. We then use a simple optical depth cri- 
terion to identify the simulation cells expected to host Lyman-limit 
absorption systems. We self-calibrate the resulting number of ab- 
sorbers to the mfp of the UVB, and present results at a given mfp 
and neutral hydrogen fraction. Our approach is fairly robust to un- 
certainties such as missing subgrid structure. 

We find that absorbers are strongly correlated with the density 
field on small scales. However, on large-scales the absorbers in- 
stead are strongly anti-correlated with the UVB, once HII bubbles 
grow enough to allow absorbers to populate these modes. There- 
fore we caution that using only a density criterion when assign- 
ing LLSs to cosmological simulations in fact neglects large-scale 
modes sourced by UVB fluctuations. 

Our ionization field algorithm assumes that a bsorbers do not 
have a large impact on reionization morphology (M cOuin n et al.l 
2007a). As such, reionization in our simulations proceeds in an 
"inside-out" fashion on large scales. Early in reionization, when the 
HII bubbles are small, the density field regulates the locations of ab- 
sorbers. As HII bubbles expand into the large-scale underdensities, 
the absorbers tend to preferentially populate the overdensities (i.e. 
filaments) of the recently-ionized IGM. Therefore the absorbers ef- 
fectively huddle around the edges of the HII regions. 

If the absorbers are indeed preferentially found on the out- 
skirts of HII bubbles, the reionization could effectively stall in the 
late stages. Reionization of the remaining HI regions might have to 
wait until the absorbers get photo-evaporated by the relatively faint 
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UVB at the edges of the HII regions (e.g. iFurlanetto & Mesingej 
2009). Alternately, new ionizing sources could reionize the remain- 
ing HI regions fro m the inside. De pending on the nature of the ab- 
sorption systems jlliev et alj r2005b1. the time-scales of both pro- 
cesses could be fairly long: of order Az ~ 1. 

If high-z absorbers are significantly neutral, they may also 
dominate the small-scale (k > 1 Mpc -1 ) 21-cm power during and 
after reionization. They smooth the contrast on moderate scales, 
and once the HII regions surpass the mfp, the absorbers add to the 
large-scale 21-cm power. Since our procedure assumes all of the 
gas inside a simulation cell is neutral, our quantitative estimates 
of the absorber contribution to the 21-cm signal can be viewed as 
upper limits. However, the qualitative trends should be robust. 

Our results should prove useful in interpreting future observa- 
tions of the reionization epoch. We also plan to extend our analysis 
to Hell reionization, where the sizable inhomog eneity of the UVB 
dPixon & Furlanettol200^ : lMcOuinn et al.l2009l) will have an even 
more pronounced impact on Hell Lyman-limit absorption systems. 
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